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Abstract —The behaviour of a high dimensional stochastic 
system described by a Chemical Master Equation (CME) depends 
on many parameters, rendering explicit simulation an inefficient 
method for exploring the properties of such models. Capturing 
their behaviour by low-dimensional models makes analysis of 
system behaviour tractable. In this paper, we present low di¬ 
mensional models for the noise-induced excitable dynamics in 
Bacillus subtilis , whereby a key protein ComK, which drives 
a complex chain of reactions leading to bacterial competence, 
gets expressed rapidly in large quantities (competent state) 
before subsiding to low levels of expression (vegetative state). 
These rapid reactions suggest the application of an adiabatic 
approximation of the dynamics of the regulatory model that, 
however, lead to competence durations that are incorrect by a 
factor of 2. We apply a modified version of an iterative functional 
procedure that faithfully approximates the time-course of the 
trajectories in terms of a 2-dimensional model involving proteins 
ComK and ComS. Furthermore, in order to describe the bimodal 
bivariate marginal probability distribution obtained from the 
Gillespie simulations of the CME, we introduce a tunable mul¬ 
tiplicative noise term in a 2-dimensional Langevin model whose 
stationary state is described by the time-independent solution of 
the corresponding Fokker-Planck equation. 

Index Terms —Gene regulatory networks, bacteria Bacillus 
subtilis , low dimensional approximation, Fokker-Planck. 


behaviour, whose noise-induced activation yields bimodal dis¬ 
tributions. An alternative model has been proposed by Suel et 
ai. mm which includes a slower, negative influence on ComK 
levels via the expression of the comS gene, shown in Fig. [I] 
This leads to an excitable system m, whose high expression 
(competent) state is not stable, but undergoes slow decay back 
to the low-expression (vegetative) state. The long-lived state 
thus explains the bimodality and removes the need to invoke a 
mechanism for the loss of competence that the bistable model 
requires. This is the model that has been well studied and has 
received considerable attention. 
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I. Introduction 

The speed and quantity of production of functional pro¬ 
teins to respond to the needs of an organism is regulated 
by transcription factor proteins that modulate the efficiency 
of mRNA transcript generation from which the proteins are 
translated, proteases that regulate their degradation, and a 
host of other molecular agents. Using fluorescent markers in 
cell-sorting ci e m and time-lapse recording technologies 
t4l . dynamic traces of the cellular behaviour have become 
accessible, making it possible to develop a quantitative un¬ 
derstanding of the dynamics of gene regulation in single cells. 
In starved conditions, some Bacillus subtilis cells undergo a 
transformation called competence, ingesting DNA from its 
environment. Single cell fluorescent images reveal bimodal 
cell populations of ComK, a protein that is used to track 
the competent state. Bimodal probability distributions can 
be generated by an underlying stochastic dynamical system 
whose deterministic state is bistable. Maamar and Dubnau 0 
demonstrate that the positive feedback provided by ComK 
proteins activating its own transcription can generate switching 


Fig. 1. Competence circuit architecture. The competence circuit includes the 
following components: two genes comK and comS corresponding to two 
proteins ComK and ComS, respectively; and promoters P CO mK and P CO mS • 
In this figure, ComK actives the expression of its own gene (auto-regulation 
feedback) and inhibits expression of ComS (negative regulation), that in turn 
interferes with degradation of ComK. The complex of MecA, ClpP/C also 
actively degrades ComK. 

The model we study m represented as a network graph 
in Fig. |T] translates into a mathematical description involving 
multiple reaction and species types in a continuous time 
Markov process, the Chemical Master Equation (CME)|[9l. 
The mean values of the random variables are described by 
a system of ordinary differential equations. While the full 
stochastic description (high-dimensional) is faithful to the 
kinds of processes believed to occur in a cell, insights into the 
dynamical behaviour of these systems are reliably extracted 
primarily from their reduced models (low-dimensional). The 
excitable behaviour of the competence circuit 0, and its 
extended design synthesized in no) . are instances of the 
application of insights gained from low-dimensional non-linear 
dynamical systems. In mu), a two-dimensional reduction 
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of the detailed model formed the basis of partitioning the 
phase space of the model into distinct dynamical regimes, 
to which the competent and vegetative states were allocated 
places. By altering the network properties away from its wild- 
type state, the authors provided evidence of the validity of 
their modelling framework when the engineered strain of B. 
subtilis followed the model prediction, for example, entering 
an oscillatory phased. However, the two-dimensional model 
proposed in m fails to reproduce two salient features of the 
stochastic model described by the CME. The duration of the 
competence state in the two-dimensional model is about four 
hours, while those in the stochastic simulations which match 
the observations from time-lapse data last around 10 hours. 
Further, the location of the second (non-steady state) mode of 
the probability distribution of the ComK and ComS proteins 
is at a different location in the reduced model compared to the 
stochastic simulation. 

The contribution of this paper is multifaceted but subtle. We 
use the method described in [IT, 12, 131 to obtain a more ac¬ 
curate low dimensional reduction for the existing model. This 
introduces a better approximation for slow-varying variables 
involving the binding of ComS to MecA, which corrects for 
non-adiabaticity in ComS dynamics. Although the underlying 
model parameters have enough uncertainty that this change in 
the competence behaviour could also be obtained by changing 
the values of the model parameters, it is nevertheless important 
in mathematical modelling to obtain a good match in the 
observed behaviour of the chemical master equation and its 
low dimensional model. A second contribution of the paper 
is to build a low dimensional stochastic model of the system. 
This is essential to capture the competence behaviour which is 
stochastically driven. We will show that the usual assumption 
of a Poissonian noise (where the variance is proportional to 
the number of reactions) provides a very poor description of 
the actual system. In contrast, we show that using a model 
in which the variance in the noise is proportional to the 
square of the molecular number can give a more accurate low 
dimensional approximation. In particular, by fitting the level 
of the noise we show that we can obtain a Fokker-Planck 
description which closely matches that observed in the full 
CME. 

In order to build a two-dimensional reduced model, we 
identify one fast species of the two bound complexes of the 
protease MecA to ComK and ComS, the principal variables 
in the model. This is in contrast to in m, where both these 
complexes are assumed incorrectly, to have fast-decaying tran¬ 
sients. However, as the ComK is not governed by exclusively 
fast reaction, we cannot use the adiabatic approximation [14). 
Thus, we use an iterative scheme proposed in [Qj] 112, [13] 
to describe the low-dimensional slow manifold as a function 
of the ComK and ComS variables. The reason for this is 
because the dynamics of the whole system breaks up into 
different time scales. Thus, although the individual species 
do not have very different time scales, there still exists 
one slow time scale describing the trajectory of the system 
after the transition to competence. This trajectory is uniquely 
characterized by the values of ComK and ComS. By doing 
this, we find that this model captures the dynamics of the 


decay of the competent state accurately. In order to account for 
the noise-driven transition to the competent state, the standard 
assumption is that for each reaction that contributes to the rate 
of change of a dynamical variable, the equality of the variance 
in the number of reactions with its mean, that is true for the 
Poisson process updates, can be used to replace the Poisson 
processes by independent Gaussian noise contributions with 
the same mean and variance. This assumption underlies the 
Eangevin approximation ElEUE) to the CME, and we find 
that it fails to reproduce the distributions obtained from the 
Gillespie simulation. Therefore, we introduce a noise term in 
the Langevin description that reflects the ratio of variances of 
the ComK and ComS distributions at the steady state. We then 
tune the magnitude of the noise so that the stationary distribu¬ 
tion, as computed from the solution of the time-independent 
Fokker-Planck equation, gives rise to a bimodal distribution of 
the ComK-ComS variables that is qualitatively similar to the 
marginal distribution computed from the Gillespie simulation. 


II. A MODEL FOR EXCITABLE DYNAMICS FOR 
EXPRESSION OF ComK 


In this section, we describe the chemical reactions provided 
in n m that give rise to the behaviour sought to match 
the observation of stochastically triggered transient compe¬ 
tence events that are identified with high expression levels of 
ComK. This set of chemical reactions describe a continuous 
time Markov process, called a Chemical Master Equation 
(CME) 0 which updates the probabilities of the number of 
molecules of each species in the system. The averages of these 
molecular numbers over multiple realizations is described by a 
set of ordinary differential equations (ODEs), called reaction- 
rate equations (RRE) El whose phase portrait best illustrates 
the excitable dynamics mm. 


A. RELATED WORK: The reactions in the discrete model 
(CME) 

The competence circuit includes the following components: 
The two principal proteins modelled are ComK and ComS and 
their corresponding promoters are P CO mK and P CO mS • ComK 
inhibits expression of ComS (rate g) that in turn interferes 
with degradation of ComK by competitively binding to MecA 
(rate & 13 ); the complex of MecA also actively degrades 
ComK (fen). The transcription of mRNAs occur in a basal 
(rates ki,k±) and regulated fashion (rates f, g ), the stochastic 
reactions are described as follows (for simplicity, we denote 
the species concentrations [mRNA corn K\, [mRNA corn s ], 
[ComK], [ComS], [MecA], [MecAx], [MecAs ], as Rk , 
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PcomK an d PcomK are constitutive and regulated promoters of 
ComK, respectively. Rk and Rs are mRNA molecules from 
which proteins ComK and ComS are translated, respectively. 
The symbols above the arrows denote the probabilities of 
reactions in unit time. The gene regulation functions /, g are 
of the Hill type: 
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The rates in the above equations are defined via the average 
rates of change as they appear in the RRE, and are functions 
of concentrations — the number of molecules per unit volume. 
For bimolecular reactions, the probability of one molecule to 
find another is inversely proportional to the cell volume and is 
tracked by the parameter Q. By using the same chemical con¬ 
vention used in then Q « InM (nano-molar); therefore, 
we can treat the concentrations of species in the same way as 
their molecular number by measuring concentrations in units 
of nM. 

The deterministic description of the system are described 
by the following differential equations (the discrete model 
parameters are given in Table [I]): 
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Fig. 2. Trajectories generated by the Gillespie algorithm are shown in thin 
lines. Also shown are the nullclines of ComK (dash line), ComS (dotted line), 
and fixed points which are the intersections of nullclines. Stable fixed point is 
denoted by full circle, saddle point by empty circle, and other unstable fixed 
point by empty square. 


The trajectories generated by the CME can be simulated 
by the Gillespie algorithm m. Fig. [2] shows such trajectories 
(after the initial transient) in thin lines plotted on the log-scale 
phase plane, for the model parameters and initial conditions 
in Tables HQ In addition, we plot the ComK and ComS 
nullclines defined by dK/dt = 0 and dS/dt = 0 with dash 
and dotted lines, respectively. Their intersections are the fixed 
points of the system dynamics. There are three fixed points, 
one of which is stable (full black circle), corresponding to low 
levels of K\ the middle one is a saddle fixed point (empty 
circle). The intersection of nullclines at the highest value of 
K corresponds to an unstable fixed point (empty square) and 
explains why the trajectories do not appear in the centre of the 
figure. The dense area around the stable fixed point is called 
the vegetative state and the trajectories that extrude into large 
ComK expression levels correspond to the competent state of 
the organism. 

There are two critical properties of the system dynamics we 
are interested in: the first one is the competence duration which 
is the time spent by a trajectory in the state where ComK 
has high level of expression, which we take to be ComK 
numbers above 10 4 (this threshold value is the same as that 
was used in ( 6 ); therefore, our results are directly comparable 
with their findings); the second is the stationary probability 
distribution which describes the probability of the system 
being at a particular state. In order to understand the system 
dynamics, a phase plane analysis was carried out in (Slz) as 
is customary 03 . However, the reduced 2D model, based on 
an adiabatic approximation where all variables except K and S 
proteins are eliminated, has a time course in the competence 
regime that is far from that produced by the full model in 
( II.2| ). This is shown in Fig. [ 6 ] where the competence duration 
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TABLE I 

Parameters of the discrete model (source from ( 3 ) 


ki 

0.00021875s- 1 

k 7 

0.005s- 1 

k\2 

0.05s- 1 

k 2 

0.1875s- 1 

kg 

10- 4 s- 1 

k\ 3 

4.5 x 10- 6 niW- 1 s- 1 

k 3 

0.2s- 1 


0.005s- 1 

k- 13 

5 x 10- 5 *- 1 


Os -1 

kio 

10 4 s 1 

k\A 

4 x 10“ 5 5- 1 

k 5 

0.0015s- 1 

kn 

2.02 x W-^nM-'s- 1 

kk 

5000nM 

k 3 

0.2s- 1 

k- 11 

5 x 10 _4 s _1 

k s 

833nM 


TABLE II 

Initial conditions 


r pconst 1 

E comK -1 

InM 

R k 

1000 nM 

M s 

100 nM 

r pconst] 
\-*comS J 

InM 

Rs 

1000 nM 

K 

1000 nM 

\Pcom K.\ 

InM 

A 

300 nM 

S 

100 nM 

\PcomS ] 

InM 

M k 

100 nM 




is a factor of two smaller in this adiabatic model than in the 
full 7D model. In the next section, we present a method of 
approximating our full system by a two-dimensional one while 
retaining the system’s competence duration. In this method, 
we first eliminate variables which are seen as fast variables in 
order to reduce the 7D system down to a 3D system. Next, we 
continue to reduce the 3D system to a 2D system by using an 
iterative procedure to capture the slow manifold of the system. 


B. THREE-DIMENSIONAL APPROXIMATE SYSTEM 


In order to do the model reduction, we first realize that the 
dynamics of mRNAs are much faster than that in the protein; 
therefore, we can eliminate the mRNAs by setting them to their 
steady state values (dRK,s/dt = 0). We restrict our attention, 
in this section, to the competent regime. Furthermore, since the 
decay rate fa^ of Mk is 500 times and 10 times faster than the 
proteins and mRNAs respectively (notice that even though the 
binding/unbinding rates of MecA to ComK ( fai/k-u ) and 
ComS (& 13 A- 13 ) are similar, the decay rate of the complex 
Mk is much faster than that of Ms, thus the dynamics of 
binding of MecA to ComK is treated different from its binding 
to ComS). We therefore can approximate the dynamics of Mk 
by slaving it to the other variables. That is, we replace it by the 
steady state value Mk * obtained from solving « 0, 

which gives Mx* = AK/T with = fe ~ 1 ^ fcl2 . Using the 
conservation equation Mt = A + Mk + Ms, we obtain (for 
Q := A + Mk) 


M k 


K 

f k + K 


Q 


(There is some freedom in choosing which variables to elim¬ 
inate and which to keep. Although this leads to the same 3D 
approximation, in making a further reduction to a 2D, we find 
that using Q rather than A leads to a series of approximation 
which converges faster to the behaviour of the 3D system. 
Intuitively, variable Q represents an effective concentration of 
MecA for the ComK dynamics, the fraction that is not taken 
up by ComS). 


After making these approximations, the 3D system can be 
described by the following differential equations: 


fa 


dK 
dt 
dS _ 
dt 1 
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UK 2 


k 7 V" x ' k 2 k + K 2 
k 5 k 6 /k 9 


ki 2 KQ 

U + i? 


-hK 
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+ k-is(MT — Q) 


kinS — kisTk- 


SQ 
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5 = (*14 + *-13)(M t -Q)~ 
at 1 k + A 


(H.3) 


From now on, we will take the initial condition to be 
K = 1099, S = 564, A = 16, M K = 3 , M s = 481, 
Rk = 1, Rs = 0 (these values come from a CME simulation). 
Fig. [3] shows that the trajectories from this 3D model ( II.3 ) 
closely follow those from the 7D model l |lL2| >. We now use 
the 3D model to find an approximation using a 2D model. 
The adiabatic approximation used in 11 eo is one such 2D 
model which gave the same fixed points as the 7D system by 
construction, but which gave rise to competent behaviour that 
was smaller by a factor of 2. This discrepancy is also shown 
in Fig. [3] This implies that the approximation dQ/dt ~ 0 used 
there assumes that the species Q changes faster than it really 
does in the 7D case. Moreover, the complexes Mk and Ms are 
treated in identical ways in mu even though their decay rates 
of those species are significantly different. Also in Fig. [3j the 
7D model trajectories are close to the ComK-nullcline of the 
2D adiabatic model, yet the 2D adiabatic model trajectories are 
not. The 3D approximate model shows that these trajectories 
and the ComK-nullcline are close to each other. This means 
the 3D approximate model has captured the slow manifold 
in the competence regime. In the following section, we first 
identify the existence of fast modes in the system which 
justifies the search for a low-dimensional approximation, and 
then introduce a procedure to find one. 
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Fig. 4. The spectra of eigenvalues on a 10-based logarithm scale. The eigenvalues are computed along the average trajectory [(a)] the position of the eigenvalue 
point is defined by the angle formed by the data point at which the eigenvalue is evaluated and the vertical axis of th e polar coordinates, and the distance 
from that point to the origin. This distance is computed by taking the logarithm of the inverse absolute eigenvalue |(b)| 



Fig. 3. Trajectories in the approximate and full systems together with the 
nullclines of the 2D adiabatic approximate model. We compare the t rajec torie s 
generated by the 3D approximate and the full models given by ftl.3) and \ll.2) , 
respectively. The numerical initial condition for the integration is K = 1099, 
S = 564, A = 16, M k = 3, M s = 481, R K = 1, R s = 0. The 
trajectory generated by the 2D adiabatic approximate model is also shown for 
comparison, where it does not match that in the full system. The trajectories 
of the 3D approximate model are close to the ComK-nullcline, this implies 
that this model has captured the slow manifold in the competence regime. 


C. A TWO-DIMENSIONAL APPROXIMATE SYSTEM 


As can be seen in ( |IL3| ), there are no obvious fast variables; 
not the Q-variable — the factor (£q 4 + fc_i 3 )/fci 3 in the Q- 
evolution equation appears also in the S evolution equation, 
but its magnitude is much smaller than that of k \2 that appears 
in the K equation — and not any of the others. This suggests 
the possibility of fast modes that decay rapidly leaving the 
long-time dynamics dependent on a reduced subset, which 
we choose to be K and S. In order to verify this, we plot 
the magnitudes of the eigenvalues of the Jacobian along the 
average trajectory to see if they are widely separated. These are 
plotted in polar coordinates in the K-S plane in Fig. [4] where 
we have chosen an arbitrary point in the competence region 
as centre. For a dynamical system dX(t)/dt = /(X(t)) with 
variables X(£) = the linearisation of the 

evolution equations around each point X(£): d(S'K(t))/dt = 
/(X + SX.) — /(X) = J5X(t) defines the Jacobian matrix J 
with matrix elements Jij = • The Jacobian matrix for 

the approximate system is defined as, 


(K) (S) (Q) 

(X\ / d/i d/i 

dK as dQ 

T — ( d/2 d/2 d/2 

J — dK dS dQ 
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+ 13 ( M t — Q ) 

h = (^14 + fc_i 3 )(Mt — Q) — ^isr*;- 


SQ 


■K 

We can express 5X(t) using the eigenvectors of J as a basis: 
<SX(£) = J^QZ^e Ait , where ^ are the right eigenvectors 
corresponding to eigenvalues A i9 and c % are the components of 
X(t) along Vi. Negative eigenvalues with large absolute values 
imply that deviations decay rapidly, and a widely separated 
set of eigenvalues enables us to eliminate these fast decaying 
modes. Fig. [4] shows a plot of eigenvalues computed along the 
average trajectory for the 3D system described in S on a 
10 -base logarithm scaled polar coordinates. 

In order to compute the average trajectory, we sample 
the data from the Gillespie simulation of the full system; 
we then divide the trajectories into bins defined by their 
angular position with respect to an origin chosen in the centre 
of the trajectories (Fig. |4(a) }. Next, we compute the mean 
eigenvalues for the Jacobian matrix computed in each bin as 
shown in Fig. |4(b)| In this figure, a particular eigenvalue is 
plotted in such a way that the distance from it to the origin is 
computed by taking a 10 -base logarithm of its inverse absolute 
value. 

It is clear that the 3 eigenvalues are separated from each 
other during the excitable state back to the vegetative state, 
making possible a reduction to a lower-dimensional system. 
In our case, the most negative eigenvalues are about 10 times 
as large as the others in absolute value, implying the existence 
of a low-dimensional attracting manifold. However, there also 
exists positive eigenvalues marked in Fig. [4] which is the 
hallmark of an excitable system. For this reason, the whole 
space is divided into regions which are defined by positive 
and negative eigenvalues; the regions where the positive eigen¬ 
values are found are demarcated by angles a , /3. In these 
regions, the trajectories could diverge along the direction of 
the eigenvectors, making the time series of the reduced and 
high-dimensional model different from each other, while the 
phase plane trajectories get pulled back to coincide on the 
low-dimensional manifold. 

Since we are interested in the dynamical behaviour of the 
reduced model on (AT, S ) space, we assume that Q can be 
described as a function of K and S, that is we assume Q = 
Q(K, S ). This means the dynamics of the system always lies 
close to a 2 D manifold in (K,S,Q) space and its velocity 
is uniquely determined by K and S alone. As a r esult, we 


have + ifQijt- Plugging this back into dll.3b we 


obtain Q = F(K , (The functional form of F can 

be found in Appendix [Bj. 

In order to estimate function Q , we use an iterative pro¬ 
cedure CEUElCCa to define a sequence of approximations 
Q n , starting with a trial function Qq(K,S) = 0 . Iteratively, 
we compute Q n +\ = n = 0 , 1 ,.... 

Each step of the iteration Q n can be put back into pi. 3) 


to obtain a 2D deterministic approximate model. Numerical 
experiments show that Q n converges rapidly and even Q 2 
gives a very good approximation. In particular, Fig. [5] shows 
the three different 2D models corresponding to the three 
first approximate functions of Q (Qi,Q 2 ,Qs) compared with 


the 3D model. It turns out that Qi is the same expression 
as that obtained by setting dQ/dt « 0, the 2D adiabatic 
approximation. It is clear that the third approximation Q% very 
closely fits the 3D approximate model; therefore, we will take 
Qs as the deterministic approximation to the full system. 



Fig. 5. Trajectories generated by the 3D model and its 2D approximations 
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Fig. 6. Competence duration using Q 3 (the 2D approximation) and Q 1 (the 
adiabatic approximation) in comparison with the 3D approximate model. Note, 
the curves of Q 3 and the 3D approximation fall on top of each other. 

Fig. [ 6 ] shows a comparison of the competence durations 
between the naive adiabatic approximation introduced in Elia 
and the iteratively produced model where Qs{K^S) is used 
in the AT, S evolution equations. Evidently, the competence 
duration in the 2D approximation (Q 3 ) is about ten hours 
which agrees with that in the full system whereas this duration 
is only roughly four hours in the adiabatic approximation. 
This significant discrepancy implies that the naive adiabatic 
model provides a poor approximation of the original system. 
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However, even though the 2D approximate model gives a much 
better approximation where the competence duration has been 
preserved, we still need to verify if this model gives similar 
stationary probability distribution by looking for a reduced 
stochastic model, which we address in the following section. 


III. THE REDUCED STOCHASTIC MODEL 


The CME was set up to capture two key properties that can 
be measured experimentally using fluorescent tagging. These 
are the time spent in the competence state and the probability 
of becoming competent. For our reduced 2D approximation to 
provide a good model of the full system, we want it to capture 
these same properties. To accomplish this, we want to make 
the probability density function (PDF) in the K-S plane to be 
similar for both the 2D approximation and the CME. Having 
obtained a procedure for obtaining a 2D model for the time 
duration, we now move on to construct a set of stochastic 
differential equations (SDE), also called Langevin equations 
0E), that reproduce the histograms of the AT, S proteins 
obtained from the Gillespie simulation. For the duration of 
competence events, we assumed that our system — the full 
model and the approximation — was initialized to start in 
the competent regime. The introduction of stochasticity is 
designed to drive the transition from the stationary vegetative 
state to the competent state via fluctuations in species numbers, 
since the process of transcription and translation are stochastic 
events m. 

The standard way of reducing a CME to a Langevin de¬ 
scription can be described in two steps 0. Firstly, we replace 
the Poisson distributions implicit in the CME description with 
Gaussian distributions for the multiplicative fluctuations in the 
Langevin equation. Secondly, we adiabatically eliminate the 
fast reactions (as before) to obtain a reduced set of equations 
describing the dynamics of slowly varying species. For the 
competence regime, the elimination of the mRNA species 
did not introduce any errors in the dynamics. However, we 
find the variance in the noise is not linear in K\ therefore, 
the standard assumption that the variance is proportional to 
the number of reactions is wrong. Indeed, Fig. [7] shows the 
significant discrepancy in the dynamical behaviour between 
the 2D Langevin model and the full system near the steady 
state. The reason for this, we believe, is: the transition to 
the competence state is driven by small number of mRNA 
molecules subject to the K 2 non-linearity in the activating 
dynamics (see the Rk evolution in ( |II.2| ), for K kj . c ), 
rendering the replacement of the Poisson by a Gaussian 
distribution suspect. Furthermore, the tails of the distribution 
of Rk molecules are drawn away, not back to the attracting 
steady state, via the repulsive dynamics that is characteristic 
of the excitable system (see the range of positive eigenvalues 


in Fig. 4(b)). In addition, even at the steady state fixed point 
the Jacobian is a non-normal matrix, and its non-orthogonal 
eigenvectors give rise to large transients driving transitions to 
the competent regime in this excitable system. Consequently, 
the elimination of the mRNA species of ComK leads to loss 
of the fluctuations needed to drive the competence behaviour. 
This is further reinforced by the observation that the linear 


noise (Gaussian) approximation around the steady state fails 
to account for the tails of the distribution around the steady 
state. 



K 

Fig. 7. Comparison between the full system and the 2D Langevin model 
in terms of probability density function (PDF) of ComK. In both models, 
we sample the simulation data for K and S satisfying 0 < K < 200 and 
0 < S < 1000. The PDF of ComK for each model is then computed and 
compared. 


By tracking the noise-driven transition in a 2-species feed¬ 
back circuit which exhibits the similar behaviour around the 
steady state, we found that the variance of the noise term 
is proportional to the square of the molecular number (see 
Appendix [C]). This finding therefore motivates us to introduce 
a Langevin model with a multiplicative noise extension of the 
2D approximate ODEs from the previous section. This 2D 
stochastic model can be described as follows, 


dK = f k (K, S, Q 3 (K, S))dt + ii k dw k 
dS = f s (K, S , Q 3 (K, S))dt + n s dw s 


(HI.l) 


where the functional forms /& and f s are given in (II.3), we 
introduce Wiener processes dwk and dw s , and set fik = cr^K, 
fi s = cr s S. This is because we are only interested in the 
tail of the stationary probability distribution where the noise 
terms are supposed to be proportional to K and S. The 
magnitudes cr/c, a s of the noise terms are chosen to obtain the 
dynamical behaviour qualitatively similar to that of the CME. 
Additionally, the initiation probability of competent events 
computed from the stochastic model should also quantitatively 
be preserved. However, this probability is very sensitive to 
the switching behaviour driven by the noise terms. In other 
words, a slight change in coefficients cr/e, cr s leads to a 
significant change in the initiation probability. This is because 
of the exponential sensitivity of the tail of the probability 
distribution to the noise-driven switching state in our particular 
circuit ED. Consequently, we have tried the simulation with 
different values of crcr s and found that the stochastic 
model which provide a relatively good approximation has 
cik = 0.008, (j s = 0.005 (see Appendix We show in 
Fig. [8] trajectories generated by the Euler-Maruyama method 
1201 from the 2D stochastic model described by ( |III.1| ), which 
are similar to those generated from the CME by the Gillespie 
algorithm (Fig. [2]). In order to obtain the stationary probability 













distribution, we solve the Fokker-Planck equation m which 
is described as follows: 

= ~^f k (K,S,Q 3 (K,S))P(K,S,t) (HI.2) 
--^f s (K,S,Q 3 (K,S))P(K,S,t) 

+ \ S ’ t] + \^ 2 e P ^ *) 
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Fig. 8. Trajectories generated by the 2D stochastic model given by with 
coefficients cr & = 0.008, cr s = 0.005. The nullclines of the 2D approximate 
model ( Q 3 ) are also plotted. 
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dimensional deterministic dynamical system, and was used as 
the basis of the stability analysis in (7). We find, however, that 
the fast reactions do not segregate the fast and slow species 
in the model description. We propose a suitable functional 
dependence between the variables that allows the system to be 
described by a slow set of variables that adequately matches 
the time scales of the competence regime, unlike the standard 
adiabatic approximation of Q. 

In accounting for the noise-driven transition to the com¬ 
petence regime, we needed to put back the noise that the 
mRNA species contributes to this process. At the vegetative 
steady state, the balance between the protein production terms 
and the linear decay term might suggest that the variance 
of the Poisson distribution be proportional to the protein 
concentration. By explicit simulation of the ComK mRNA- 
protein subset, we find this does not hold, and instead find that 
the standard deviation is linear in the protein concentration. 
Using this as a guide, we introduce a 2-dimensional stochastic 
differential equation with a Wiener process proportional to 
a linear term in the concentration, and fit the constant of 
proportionality by minimising the distance between the distri¬ 
butions obtained from the chemical master equation and the 2- 
dimensional approximation. This approximate 2-dimensional 
system could then form the basis of future studies on the 
dependence of first passage time distributions on the system 
parameters, and other theoretical characterisations. We propose 
that this indirect way of approximating the behaviour of a 
high dimensional chemical system by a low-dimensional set 
of stochastic differential equations could be used as a generic 
method when direct ways of model reduction fail. 


Fig. [9] shows the probability density function (PDF) of the 
2D reduced model in comparison with that computed from 
the full model. As we can see, both models produce similar 
bimodal distributions that are characteristic of the cell counts 
in the vegetative and competence states obtained in (6). 

IV. CONCLUSION 

In this paper we have systematically conducted a series of 
customary approximations of the chemical master equations of 
the excitable model of competence and found many to be not 
faithful to the dynamics of the full model. First, to obtain 
the ODE description of (3 we replaced the stochastically 
varying quantities by their mean values and further assumed 
that the averages of nonlinear functions (f(x)) of a stochastic 
variable x can be replaced by non-linear functions of the 
averaged variable f((x}). The reduction of the ODEs thus 
obtained to a two-dimensional set is usually accomplished 
by dividing the dynamical processes into a set of fast and 
a set of slow reactions. If the variables are also divided in a 
similar way into a fast and slow set, a fast-slow decomposition 
of the dynamics is easily facilitated. This allows the slow 
variables to be viewed as unchanging when considering the 
fast dynamics; the slow variables are only affected by time 
averages of the fast variables. This adiabatic approximation 
— replacing fast variables by their equilibrium values — turns 
out to be useful in retaining the same fixed points as the high 


Appendix 

A. Derivation of Eqs. S 

The first model reduction step would be eliminating both 
mRNAs in Eqs. dIT2l > by setting dRx,s/dt = 0, we obtain 


Rk = 
Rs = 


k\ + 


kp+K* 


k 7 

k 3 


1 + (. K/k s ) 5 

On the other hand, we also have, 

K -M 


kg 


Mk = r„ + k 

Mg = Mt — Mk — A ~ Mj 1 — Q 

K 


A = Q — Mk = Q — 


r k + K 


Q = 


K 


Q 


From Eqs. ( |A.1| ), ( |A.3| ) and ( |A.5| ), we have 

dK 


—j— = —kuKA + k-nMx + ksRx — kgK 
dt 


= -kuT 


ii-L k j 


KQ 


+ k- 


h 

k 7 
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f K 
k 2 K 2 
k 2 k + K 2 
k 2 K 2 
ki + K 2 


it 


KQ 
r k + K 

-k 8 K 

k\ 2 KQ 

f V+K 


(A.l) 

(A.2) 

(A.3) 
(A .4) 
(A.5) 

(A.6) 
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Fig. 9. Probability density functions (PDFs) of the phenomenological stochastic 2D model [(a)] and the full discrete model [(b)] 


where fci 2 = fc_n — 
CO), we have 


knT k . Similarly, from (A.2), (|A.4|) and 


dS 

dt 


—ki^SA + + ^6^5 — fcioS 


k 5 k 6 /k 9 
1 + (iT/fc s ) 5 


— feioS 1 — fci 3 r fe 


SQ 

r k + K 


+ A:_i3(Mt — Q) 


(A.7) 


(A.8) 
(A.9) 


Finally, 


dQ 

dt 


dM s 

dt 


(A. 10) 


(*- 13 + ku)Ms — ki^SA 

(*-i 3 + * 14 ) (m t - 0 ) - *i 3 r fe — 

f /c + A 


From ( |A.6| ), ( |A.7| ) and ( |A.1Q| ), we obtain B- 


Iterative Procedure For Finding Q 

In order to find the explicit form of Q , we rewrite the 
equation as follows: 


(&14 + k-is)(M T — Q) — fci 3 rfc 


SQ 


dQ (k% . 

= dK U ( 1 + 


* 2 ir 2 

fc2 + /f 2 


r fe + K 

ki 2 KQ 
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— ksK 


dQ ( k 5 k 6 /k 9 SQ 

+ dS )IT(W " 10 “ 13 ‘ ivn? 


+ k-iz(M T — Q ) 


hence, 


Q = F 


dQ dQ 
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C. Reduction of a 2-Species Model to Track Noise-driven 
Transition 

As we mentioned early in section |In| the usual Poissonian 
noise in the Langevin approximation turned out to be incorrect. 
In order to address the source of this issue, we will focus our 
attention on the behaviour of the system near the stable fixed 
point and the transition beyond the intermediate unstable fixed 
point to the competent state. To do this, we will be looking at a 
much simpler noise-driven switching circuit which is extracted 
from the original system. This is done by just looking at the 
dynamics of two variables ComK and mRNA while ignoring 
the effect of the other variables, this will give us a bistable 
model. Even though the behaviour of this bistable model is 
different from the full model in a long period of time, the 
dynamics of the system near the fixed point should be similar 
in both models. Hence, if we can capture the noise near the 
transition which is believed to account for the tails of the 
distribution around the steady state, we then can construct a 
stochastic model which describes the similar behaviour as that 
in the full system. 

The two-species model is obtained by only considering 
the chemical reactions where a protein activates its own 
transcription as follows: 
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TABLE III 

Model parameters 


ki 

0.00021875s -1 nM 

/C 2 

0.1875s- 1 

k 3 

0.2s- 1 


0.005s- 1 

k 5 

3.2 x 10- 5 s- 1 


1.4704 x 10 -4 s -i 

kk 

5000 nM 


PSmK^PS + mRNA comK 

PcomK WCornKlM^ + m R N A comK 

mRNA comK mRNA comK + ComK (A. 11) 

mRNA comK ^> 0 

ComK-^ 0 


where f ([ComK], k 2 ,kk) = ' 

The first two reactions represent how much mRNA corn K 
is produced from the binding of protein to the promoters on 
DNA. The next reaction shows how much protein ComK is 
synthesized from mRNA corn K- The fourth and fifth reactions 
represent the linear degradation of the mRNA and protein, 
respectively. In fact, this model is simplified from the 7D 
model by setting the variables MecA and MecAx to their 
steady values. Moreover, the system exhibits bistability and 
transition from a low to a high expression state of ComK 
which is driven by noise in mRNA levels. 

We denote the protein and mRNA as K and m, respectively. 
As a result, the deterministic differential equations for this 
model are described as follows: 


dX 

dt 

dm 

dt 


= &5 + k 3 m — k^K 


— ki + 


k 2 K 2 
h 2 + K 2 


k^m 


(A. 12) 


The model parameters are given in Table [III] here we introduce 
the parameter k$ (k$ = 3.24 x 10 -5 ) in order to keep the 
structure and the position of the fixed points the same as that 
in the original 7D system. 

In this model, we found that the 2D Langevin simulation 
breaks down due to the very small number of mRNA pop¬ 
ulation having been driven to negative values. On the other 
hand, we notice that the mRNA lifetimes are shorter than 
protein lifetimes = 34 >> 1), this therefore motivates 
us to adiabatically eliminate this small variable to obtain a 
ID Langevin model. However, the simulation result shows 
that this model remains around the steady state and never 
goes to the high expression regime. This means the adiabatic 
approximation does not capture the noise-driven transitions 
in this model, and the fluctuation in the mRNA which is 
generated in the 2D Langevin model significantly contributes 
to the switching behaviour of the system. Moreover, the 
mRNA reduction is also problematic due to the fact that the 
time scales of mRNA and protein are not completely separated 



Fig. 10. The square of size of fluctuation (X 2 ) is well fitted by a quadratic 
curve. This means the noise term in the Langevin model should be propor¬ 
tional to the number of ComK. The discrepancy between the variance of 
Poisson noise (tr 2 ) and the empirical noise (X 2 ) is also shown. 


(despite the decay rate of mRNA being 34 times faster than 
that of protein, the production rate of mRNA (ks = 0.2) 
is about 1000 times faster than the decay rate of protein 
(&6 = 1.4704 x 10 -4 )). Even though the mRNA reduction 
issue can prevent us from coming up with a good approximate 
model, there is a way around it. In particular, we will estimate 
the size of fluctuation near the steady state in the 2D Gillespie 
model whereby we would hope to construct the correct noise 
for the stochastic model. 

In order to estimate the size of fluctuation, we run Gillespie 
simulations using the reaction scheme given by |A.11| with the 
initial condition given by the fixed point. For each run, we stop 
the simulation as soon as the molecular number of protein 
exceeds 500. This is the threshold over which the system 
enters the high expression state. Next, we collect and put the 
simulation data into 500 separate bins according to different 
values of protein K (notice that we are only interested in 
the values of the protein and mRNA, not the time step). In 
particular, each bin i = 1,2,... 500 contains a particular value 
of protein Ki and a set of all possible values of mRNA with 
respect to TQ ( we don’t count the frequency of mRNA). Let 
Li be the total number of mRNA values in bin i, then the 
value of an instance of mRNA j belonging to bin i is denoted 
as mij where j = 1,2,... Li. For each bin i = 1,2,... 500, 
we compute the expected change in protein Ki in time step At 
denoted as ATQ that is determined by the propensity functions 
in which the protein gets involved. According to this, the 
expected change in Ki given m^ in a time interval At denoted 
as A K\ can be estimated as A Kj = (fcs + k^m^ — k^K^At. 
In our case, we take At to be the same as that in the ID 
Langevin model. Since At is the same in both Gillespie and 
ID Langevin models, the only comparable term would be 
&5 + k^mij — k^Ki. Thus, we can ignore At and re-define 
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A K\ as follows: 

A K{ =k 5 + k 3 mij - k 6 Ki (A.13) 

where j = 1, 2,... L^. The variance of A Ki is then computed 
by the following equation: 

4 = f E(A K\ - (A Ki)) 2 (A. 14) 

1 3 = 1 

Here, (A7Q) = A. A K\. On the other hand, since CME 

can be described by a Poisson process; therefore, the variance 
in protein caused by this process is given as follows: 

£i = fa + k 3 (mi) + k 6 Ki (A. 15) 

Here, the mean of mRNA for each bin i is computed as 
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(A. 16) 


As a result, the size of fluctuation for this particular data bin 
is given by 


s = a A; 2 + 


i t a A K, 


(A. 17) 


By estimating the size of fluctuation in ComK (E), we found 
that this quantity is approximately proportional to K. Indeed, 
Fig. 10 shows the estimated square of size of fluctuation (E 2 ) 


in ComK which we can fit by a polynomial fitting curve. The 
fitting is done by minimizing the sum of the squares of the 
deviations of the data from the empirical curve (least-square 
fit). Since we are only interested in fitting the part of the curve 
which account for the tail of the probability distribution, we 
will do the fitting for 100 < K < 500. In Fig.|T 0 | the empirical 
curve is well fitted by a quadratic curve which is defined as 
follows: 

Vo = b 2 K 2 + b\K + b 0 


Fig. 11. Distance as a function of cr and cr s . 


D. Parameter Optimization Procedure 

As mentioned earlier, we will try to optimize the parameters 
c J/c, cr s such that the stationary probability distribution obtained 
from solving the Fokker-Planck equation ( |III.2 } close to that 
computed from the full discrete model. To do so, we use 
the Jensen-Shannon divergence ED which is a smoothed 
version of the Kullback-Feibler divergence [22] to measure 
the distance between the two distributions, says P and Q. 
For each pair of parameters (cr/-, a s ), we compute the Jensen- 
Shannon divergence D(P, Q ) and sort them in descending 
order, we then choose the pair of parameters corresponding to 
the smallest distance. The pair of parameters chosen for our ex¬ 
periment satisfies 0.005 < < 0.02 and 0.001 < a s < 0.02. 

Fig. [TT] shows the distance between the PDFs in the full model 
and stochastic model as a function of the noise terms. In our 
experiment, we found that the stochastic model provides the 
best approximation with = 0.008, a s = 0.005. 


where 62 = 1.1 x 10 6 , 61 = —0.00014, bo = 0.15. Also in 
Fig. 10 the variance of Poisson noise (cr 2 ) when adiabatically 

fci+ k i K * 

eliminating the mRNA by setting m = - k t +K is shown 




for comparison. This quantity is estimated as follows: 


v p = fci 


h 


k 2 K z 


k 6 K 


As we can see in Fig. [T0| there is a significant difference 
between the variance of Poisson noise and the empirical noise 
(E 2 ). This implies that the Possion noise does not capture the 
correct noise in the Gillespie model which can only be fitted 
by a much bigger noise, and the magnitude of this noise is 
demonstrated by the quadratic fitting curve. We also notice that 
the empirical curve does not grow linearly to K\ therefore, it 
should not be approximated by a linear line. This means the 
variance is proportional to the square of the molecular number 
ComK. 
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